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We propose the use of two-dimensional photonic crystals with engineered defects for the generation of an 
arbitrary-profile beam from a focused input beam. The cylindrical harmonics expansion of complex-source 
beams is derived and used to compute the scattered wavefunction of a 2D photonic crystal via the multiple 
scattering method. The beam shaping problem is then solved using a genetic algorithm. We illustrate our 
procedure by generating different orders of Hermite-Gauss profiles, while maintaining reasonable losses and 
tolerance to variations in the input beam and the slab refractive index. 
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1. Introduction 

Laser beam shaping, defined as redistributing the ir- 
radiance and phase of a beam, is of great interest for 
many applications such as image processing and holog- 
raphy m , atom guiding [2j] , materials processing ^ and 
controlling random laser emission Shaping can be 
achieved using various optical apparatus, such as binary 
holograms Q , conical lenses H i] , solid state lasers Q , 
and spatial light modulators [4|. Beam shaping using 
anisotropic photonic crystals has also been reported [8|- 
[Toj . Moreover, the generation of self-healing, limited- 
diffraction Bessel-Gauss beams by 2D axicon-shaped 
photonic crystals has recently been demonstrated by 
Kurt and Turduev [ill, These promising results 

highlight the potential of photonic crystal engineering 
for the generation of beams of arbitrary profiles. How- 
ever, few solutions are available for robust integration 
of optical elements dedicated to beam shaping on pla- 
nar lightwave circuits. One of those is the use of a 
heterogeneous refractive index maps to convert a Gaus- 
sian beam to a Bessel-Gauss profile 13|. Nevertheless, 



planar-waveguide based photonic crystal slabs, consist- 
ing of air holes in a high index core, retain immense 
potential for fabrication of integrated optical elements 

The aim of this paper is to show that two-dimensional 
photonic crystals (PhC) can be engineered to achieve 
any specific beam profile required for a given applica- 
tion, while maintaining relatively low scattering losses. 
Theoretical PhC engineering involves selecting a num- 
ber of adjustable geometric parameters and performing 
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parametric optimization of a cost function related to the 
irradiance distribution of the scattered beam. Since the 
use of more adjustable parameters (usually) results in 
more diverse output profiles, a fast and accurate numer- 
ical method is needed to compute the field scattered by 
the PhC device. The speed of the method is critical since 
a large number of configurations must be tested. Con- 
sequently, resource-heavy finite-difference time-domain 
(FDTD) computations (l^ are not suited for our pur- 
pose. We rather use the typically faster multiple scatter- 
ing computations [H, [l^ . The first part of this paper is 
concerned with a description of the scattering approach. 
We present a derivation of the cylindrical harmonics ex- 
pansion of focused beams used to parametrize the wave 
incident on the PhC. This expansion is required by the 
multiple scattering formalism. 

In the latter part of this paper, we detail the pro- 
posed PhC devices and the optimization scheme used. 
Like Vukovic et al. [l^, we choose a basic photonic 
lattice configuration and allow individual scatterers to 
be present or absent as the only adjustable parameters, 
thereby enabling a binary encoding of the configura- 
tion space and the use of the standard genetic algorithm 
(GA) to find the configuration best suited to our purpose 
|2ll . |22| . Our results show that the optimization strate- 
gies presented in [13] can be advantageously used to de- 
sign an integrated beam shaping device. To illustrate 
this, we present engineered configurations allowing the 
generation of two different Hermite-Gauss beam profiles 
with great accuracy, and discuss the power conversion 
efficiency of the proposed devices. 

2. Scattering of complex-source beams by PhCs 

This section establishes the theoretical framework used 
to compute the field scattered by a finite PhC slab. A 
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generic two-dimensional PhC consists of an array of air 
holes in a planar dielectric waveguide, with a lattice con- 
stant of the order of the operation wavelength [iB, [l^l • 
Since our goal is to engineer the geometric properties of 
the PhC to achieve a given beam profile, we only con- 
sider finite-size slabs. For modeling purposes, we sup- 
pose that every cylinder (hole) is infinite along the axial 
z direction. The field scattered by the cylinder array is 
then given by the solution of the 2D Helmholtz equation 



[V2 + fc2(a,^y)]u(a;,y) =0 



(1) 



where a harmonic time dependence exp(— za;i) is as- 
sumed and k — kan{x,y), where n is the refractive in- 
dex. Both TM {u = E^) and TE {u = H^) polarized 
waves can be considered. The wavefunction outside the 
scatterers can be written as a superposition of an inci- 
dent and a scattered wave, u(x, y) = Ui{x, y) + Us{x, y). 
We then seek the scattered wavefunction Us{x, y) in the 
case where Ui{x,y) is a focused beam with a Gaussian 
shape in the paraxial zone. For this purpose, the in- 
cident wavefunction is represented by a complex-source 
beam (CSB). This solution has been proposed in order to 
extend the validity of the Gaussian beam (GB) beyond 
the paraxial zone [13, ■ 

Using the Green's function of the inhomogeneous 
Helmholtz equation for a point source located in the 
complex plane at coordinates x' — ixr and y' — Q, one 
obtains the CSB solution 



Ui{x,y) = H^^\krs) 



(2) 



where Hq^'^ is a Hankel function of the first kind. The 
complex distance is given by 

^ [{y-yr + ix^xT]'/' = [y' + {x~ ^xnr]'/' ■ (3) 

The complex point-source yields a directional field ra- 
diating away from the beam waist (x = 0). The 
CSB is continuous everywhere in the real plane except 
across the branch cut connecting the two singularities at 
{x,y) — {O^xr) and {x,y) — {0,—xr). For the purposes 
of this paper, we shall restrict our attention to scatter- 
ers located in the positive x plane, referring the reader 
to [1^ for regularization strategies in the waist plane. 
Since the H^^^ function converges rapidly to a complex 
exponential, one can readily show that, for x > 0, 



Ui{x,y) Ug{x,y)eTip{kxii + in/A) 



(4) 



where 



Ug(x,y) 



■nk{x — ixr) 



exp 
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y 



2 X — ixR 



(5) 

In other words, the CSB reduces to a GB of Rayleigh 
distance xr propagating along the x axis in the parax- 
ial zone. Moreover, since the CSB is an analytical so- 
lution of the Helmholtz equation exhibiting the cylin- 
drical symmetry characteristic of the multiple scatter- 
ing method, it is the ideal parametrization of a focused 



non-paraxial GB incident on an array of cylindrical scat- 
terers. 

2. A. Expansion of complex-source beams in cylindri- 
cal harmonics 

To compute the scattered wavefunction via the multi- 
ple scattering method, one needs to expand the incident 
field on a basis of cylindrical waves centered on each in- 
dividual scatterer. This section is dedicated to the ana- 
lytic expansion of the aforementioned CSB into cylindri- 
cal harmonics. Let On) be the cylindrical coordinate 
system local to the n*^ scatterer, whose center is located 
at (X„,y„). We seek a series expansion to rewrite the 
incident beam in the following fashion 



1— — 00 

On can rewrite eq. ([2]) as 

Ui{pn,On) = H^^\k\rn - T sn\) 



(6) 



(7) 



where r„ = On) and Vgn is the vector pointing from 
the center of the v}^ scatterer to the complex source 
point at coordinates {x,y) = {ixR,0). To uncouple r„ 
and ran m the argument of Bessel functions, we apply 
Graf's addition theorem This leads to the following 
expansion coefficients, similar to those found in [28j 



where 



and 



all = (-l)'i?i(fcr,„)e-^'^ 

Xn - iXR 
cos /i = . 



(8) 
(9) 
(10) 



For comparison, the expansion coefficients for a plane 
wave (PW) incident from the —x axis are given by 



■likX„ 



(11) 



The sole knowledge of the a^j expansion coefficients of 
the incident beam allows the use of the multiple scatter- 
ing method. In a nutshell, one writes the scattered field 
as a sum of cylindrical waves centered on each individual 
scatterer 



(12) 



n / 



The matrix equation connecting the expansion coeffi- 
cients can be written as Snia^i ~ 2^n„'^n'/'i with 

T^!n' = Snn'Sw - (1 - ^nn' )e'^'' (fci?nn' )««/ 

(13) 

where Rnn' is the center-to-center distance between scat- 
terers n and n', 4>n'n is the angular position of scatterer 
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Fig. 1. (Color online) Basic photonic lattice configuration 
{Na = 104). To generate a desired beam profile, defects can 
be present or absent. We impose a vertical mirror symme- 
try, resulting in 2^^ possible configurations. The dotted line 
indicates the plane used for the computation of the desired 
beam profile. 



n' in the frame of reference of scatterer n and Sni is a 
constant resulting from the application of electromag- 
netic boundary conditions. Further details are given in 

[ia[i3. 

Remarkably, except for the computation of a cylin- 
drical function, no supplementary numerical cost is in- 
volved in computing the scattered wavefunction in the 
case of an incident CSB rather than an incident PW. 
Indeed, the core operation of the multiple scattering 
method involves computing the 6„/ coefficients via a ma- 
trix inversion, whose computation scales as the square 
of the number of scatterers Ng , regardless of the shape 
of the incident beam. However, a simple analysis shows 
that the convergence of (|8]) is limited to a disk not in- 
tersecting or touching the branch cut between {x, y) = 
{0,Xfi) and {x,y) = {0,~Xfi). In other words, scatter- 
ers must not intersect or touch the branch cut for the 
expansion to be used in scattering computations. This 
restriction is not present in the case of an incident PW. 
It does not restrict the scope of our computations since 
we position all scatterers in the -fx half plane. 

3. Beam shaping computations 
3. A. Problem definition 

The objective is to find a PhC configuration which, when 
illuminated with a CSB, produces a scattered wavefunc- 
tion that matches a desired irradiance profile in a given 
plane. Let y) be the desired output wavefunction. 
The beam shaping problem can be formulated as the 
minimization of the following integral 



J \\u{xo,y)\'^ - \u{xo,y)\'^\dy 
J \u{xo,y)\'^dy 



(14) 



where Xq is the location of the target plane. This is 
equivalent to minimizing the root sum of squares (RSS) 
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Fig. 2. (Color online) Band structure for a square lattice of 
air holes of diameter D = 0.6A in a dielectric medium with 
refractive index 2.76. The location of the partial bandgap 
is shaded. Eigenmodes were cornputed using the MIT Pho- 
tonic BANDS software package [29l |. 



of irradiance variations at a set of points of the tar- 
get plane p|. It is worth noting that we do not take 
into account the phase of the output beam, only the 
amplitude. This increases the number of "acceptable" 
configurations in the problem space, at the cost of losing 
information about the coUimation of the output beam in 
the optimization process. Large variations in the output 
phase front may result in large output beam divergences, 
although this is not critical for applications such as ma- 
terials processing P, Moreover, since backscattering 
losses are mostly unavoidable in PhC devices, imposing 
a peak irradiance value is too severe a condition for the 
optimization algorithm. We rather seek a normalized 
irradiance profile, and evaluate backscattering losses a 
posteriori. 

The basic scatterer geometry (fig. [Ij is a variation of 
that presented in j^^, HO] , i.e. part of a square lattice of 
air holes embedded in a medium of index n = 2.76. The 
diameter of all holes is set to D = 0.6 A, where A is the 
lattice constant. The infinite counterpart of this pho- 
tonic lattice exhibits a partial photonic bandgap for both 
polarizations in the P— X direction (see fig. [5]). Although 
the strong confinement associated with a full photonic 
bandgap is exploited in the case of waveguide design 
[soj . it is not mandatory for beam shaping purposes. In- 
deed, the purpose of the finite PhC slab is not to act 
as a Bragg reflector, but rather to redistribute the inci- 
dent beam irradiance via multiple scattering. We shall 
therefore concentrate on operating wavelengths near the 
partial bandgap to ensure relatively strong scattering. 

For definiteness, we prescribe our incident beam as 
a TM-polarized CSB given by ^ with a half-width 
wq = 2. 5 A and a wavenumber ka — 1.76/ A for a 
Rayleigh distance xr = fcoWg/2 = 5.48A. Although 
the desired output beam and target plane can be ar- 
bitrary, for illustrative purposes we have chosen to gen- 
erate Hermite-Gauss beam profiles of half-width w at 
the device output, that is 



\um{xo,y)\^ = [n^ iOfexpi-e) 



(15) 
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Fig. 3. (Color online) Convergence of the standard GA used 
to find the configuration shown on fig U) The fitness value 
reached is 1/7 ~ 47.6. 

where ^ = \/2y/w and T-LmiO is a Hermite polynomial. 
The first two orders are 

(16) 

while 7^0 (0 = 1- For simplicity, we require further that 
the half- width w of the desired beam profile be identical 
to Wo- 

We use a genetic algorithm (GA) to find the config- 
uration best suited to the generation of a given beam 
profile [lO, [2l|. The problem encoding is binary, with 
each configuration being assigned a "genotype" of length 
equal to the number of available scatterer sites. For the 
purpose of demonstration, we have targeted symmetric 
beam shapes and have explicitly imposed mirror symme- 
try of the scatterers about the j/-axis. This effectively 
reduces the problem space dimension, but the method 
is equally efficient for asymmetric beam shapes. Each 
trial configuration is assigned a fitness value inversely 
proportional to /. Populations of 200 individuals are 
generated and evolution takes place until an optimum 
is reached, typically within a few thousand generations 
(see fig. ini). We use the standard GA evolutionary op- 
erators: roulette wheel sampling, mutation probability 
Prn = 0.002, uniform crossover with probability pc = 0.2 
and elitism. It is noteworthy that the computation of 
the fitness function, which implies a matrix inversion 
and field evaluation via the multiple scattering method, 
takes only a few seconds for one generation (200 config- 
urations) . 

3.B. Generation of beam profiles and tolerance of 
configurations 

In this section, we present the best configurations found 
for order 1 and 2 Hermite- Gauss beam profiles, exhibit- 
ing a zero and a maximum on the propagation axis, re- 
spectively. Results shown on figs. |4] and [5] highlight 
the possibility to generate order 1 and 2 Hermite-Gauss 
beam profiles with great accuracy (/ < 0.05) and are 
representative of a number of calculations that we have 
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Fig. 4. (Color online) Generation of order 1 Hermite-Gauss 
beam, (a) Optimized configuration and field profile {Ns = 
41). The target plane is indicated by a dashed line, (b) 
Comparison of computed irradiance along target plane and 
desired profile (arbitrary units). This design is characterized 
by / = 0.021,7; = 0.705. 



performed. For comparison, the error on the amplitude 
profile for the PhC device reported in ^ is around 10 
%, while the error of the integrated device proposed in 
|13l | is around 5 %. This shows that our designs perform 
equally well or better than recently proposed integrated 
beam shaping solutions with respect to the profile accu- 
racy. We also stress that the method used is not limited 
to a single lattice nor to a specific output beam profile. 
For example, we have obtained profiles with similar ac- 
curacy using a triangular lattice with the same refractive 
index. 

Since our primary goal is to obtain an accurate nor- 
malized profile via GA optimization, the best configu- 
rations found do not necessarily exhibit low backscat- 
tering losses. To quantify these losses, we compute the 
efficiency rj of the best designs by evaluating the ratio 
between the electromagnetic power transmitted in the 
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Fig. 5. (Color online) Generation of order 2 Hermite-Gauss 
beam, (a) Optimized configuration and field profile {Ns = 
28). The target plane is indicated by a dashed line, (b) 
Comparison of computed irradiance along target plane and 
desired profile (arbitrary units). This design is characterized 
by J = 0.044, ri = 0.785. 



target plane and the total incident power; that is 
_ I^oo S^{xo,y)dy 



(17) 



where Xm is the location of the input plane and Sx is 
the X component of the time-averaged Poynting vector 
[Toj . The computation of rj is achieved via numerical 
quadrature. As our computations show, efficiencies of 
optimized configurations typically fall between 70 % and 
80 %. These numbers are only 10-20 % smaller than pro- 
posed integrated beam shaping devices specifically tai- 
lored for high efficiencies: references 0, J3] report effi- 
ciencies of ~ 90 %. It is therefore quite rewarding that 
our final configurations not only provide a high profile 
accuracy, but also a low loss design. Of course, if a 
higher efficiency is critical to a given application, it is 
always possible to alter the fitness function of the GA 
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Fig. 6. (Color online) Tolerance of PhC lattice configurations 
to (a) variations of the Rayleigh distance of the input beam 
and (b) group refractive index of the slab. The design values 
are indicated by a dotted line. 



to optimize for efficiency as well. 

It is instructive to examine the tolerance of optimized 
PhC configurations to variations of the design parame- 
ters. In experimental situations, the Rayleigh distance 
may vary if the input beam focusing is more or less con- 
trolled. On the other hand, the slab refractive index 
may be fixed using the effective index approximation 
|31| . To assess the tolerance to variations of these two 
parameters, we have computed the RSS integral / for 
various values of xn/ A and n around the design values, 
while maintaining all others parameters fixed (fig. [B]). 
Results show that varying the value of xr/A by ±1, one 
full lattice spacing, preserves the low value of / (un- 
der 0.10), especially in the case of the order 1 Hermite- 
Gauss beam profile. The PhC configurations presented 
are also robust with respect to the parameter n. It is 
possible to draw two observations from these computa- 
tions. First, it is not necessary to run a GA search over a 
wide range of parameters to keep the fitness of the PhC 
designs within acceptable limits of performance even if 
some parameters are only approximately known in ex- 
perimental applications. Second, the results show that 
the fabrication of a PhC based integrated beam shaper 
operating in the infrared (Aq ^ 1500 nm, A ~ 500 nm) 
is well within reach of current fabrication techniques. 
Indeed, devices operating in that regime have been suc- 
cessfully fabricated in silicon-on-insulator material using 
UV lithography [H-ll^. 
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4. Conclusion 

In this paper, we have presented a general design 
method based on a genetic algorithm for beam shap- 
ing using integrated two-dimensional photonic crystals. 
Paramctrization of the incident Gaussian-like beam was 
achieved using the CSB solution of the Helmholtz equa- 
tion. The cylindrical harmonics expansion of the inci- 
dent CSB allows for the use of the multiple scattering 
method to compute the field scattered by the PhC slab. 
This method enables fast computation of the amplitude 
profile of the beam scattered by individual photonic lat- 
tice configurations. 

Using this design method, we have tailored photonic 
crystal devices for the conversion of a CSB to order 1 and 
2 Hcrmitc-Gauss beam profiles. The associated beam 
shaping error (< 5%) compares advantageously to other 
known integrated solutions. We also found that over 
70 % of the input beam power was channeled to the 
output beam. Although we have used a square lattice 
and required a Hermite-Gauss profile, different lattices 
and output beam profiles can be accommodated at will. 

We have also evaluated the sensitivity of the output 
beam to variations in the depth of focus of the input 
beam and the slab refractive index. Our results show 
that integrated amplitude beam shapers may very well 
be fabricated using current technology. 
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